{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(a)(b)(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error: 0.0025\n",
      "Error: 0.0238\n",
      "Error: 0.0088\n",
      "Error: 0.0088\n",
      "Error: 0.0013\n",
      "Error: 0.0000\n",
      "Error: 0.0000\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'error VS szie')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEICAYAAABfz4NwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VfWd//HXJytLAkgIYYdAAoq4EhbrUpaxLlPFVltxX0DaKvrrOEvb30zbGafzm3H66zKCtgouuGvpry0zVm0ruEMgKC6oQNgkrIEgECBk+/z+yAEvMSE3kOTc5f18PO4j937P95x8vslN3vme77m55u6IiIikhF2AiIjEBgWCiIgACgQREQkoEEREBFAgiIhIQIEgIiKAAkEkLpnZIDOrNLPUsGuRxKFAEGklMzvHzPabWXYT2941s5nB/Wlm9omZ7TOz7Wb2QlP7HA93/9Tds9y9ri2OJwIKBElwZpYWTVtrjuHui4Ey4MpG/UYBI4FnzOzLwP8BrnH3bOAU4PnWVS/SsRQIEnfMrJ+Z/dbMys1svZndFbHtn81svpk9aWZ7gZubacs0s1+a2Zbg9kszywyOMcHMyszse2a2DXi0iTLmATc2arsReMHddwFjgMXu/i6Au1e4+zx339fMmG42s3XBbGK9mV0XtL8XnBo6fPOgviHB/bSgX3cze9jMtprZZjP7iU4nSWspECSumFkK8N/Ae0B/YDLwXTO7KKLbFGA+0AN4qpm2fwTGA2cCZwBjgX+KOEYfoCcwGJjRRClPAOeb2aCIuq4FHg+2FwMXmdm/mNm5h8OmmTF1Be4DLglmE18CVgC4+xnBqaEs4G5gFfBOE4eZB9QCBcBZwFeA6c19TpGmKBAk3owBct39Hnevdvd1wBxgakSfxe7+e3evd/eDzbRdB9zj7jvcvRz4F+CGiGPUAz9290MRxzjC3TcBrwHXB02TgU7AC8H2N4CvA2cHbbvM7OfH+Ku9HhhlZp3dfau7r4zcaGbnAT8BLnf3vY225QGXAN919/3uvgP4RaOviUiLFAgSbwYD/czss8M34H8DeRF9NjWxX+O2fsDGiMcbg7bDyt29qoVaIk8b3QA87e41hze6+4vufhkNM40pwM008Ve7u+8Hrga+DWwNFp9PPrzdzAbSsP5wk7uvbqKOwUB6sO/hr8mDQO8W6hc5igJB4s0mYL2794i4Zbv7pRF9mvoXvo3bttDwi/SwQUHbsY7R2P8D+pvZRBpmA4831SmYlbwCLARGNdPnZXe/EOgLfELDrAcz6wz8Hvilu7/YTB2bgENAr4ivSTd3PzWKMYgcoUCQeLMU2Bss+HY2s1QzG2VmY1p5nGeAfzKzXDPrBfwIeLI1Bwj+sp9Pw6LzRncvObzNzKaY2VQzO8kajAW+DCxpfBwzyzOzy4O1hENAJXD4ctJHgE/c/T+PUcdW4E/Az8ysm5mlmNmw4EonkagpECSuBNfdX0bDYvB6YCcwF+jeykP9BCgB3gc+oGGh9ifHUdI8GmYajWcHu4HbgDXAXhrC5qfu/hRflAL8LQ0zlAoaguP2YNtU4GuNrjQ6v4lj3AhkAB8Fn3s+DbMNkaiZ3iBHRERAMwQREQkoEEREBFAgiIhIQIEgIiIAtOqffIWtV69ePmTIkLDLEBGJK8uXL9/p7rkt9YurQBgyZAglJSUtdxQRkSPMbGPLvXTKSEREAgoEEREBFAgiIhJQIIiICKBAEBGRgAJBREQABYKIiAQSPhDq6p1nl37Kix9sDbsUEZGYFlcvTDseKQZPL/2Uzw7UcOHIPNJSEz4DRUSOS8L/djQzZk4s4NOKAyx4b0vLO4iIJKmEDwSAC0fmcXKfbGYvKqWuXm8IJCLSlKQIBDPjzkmFrCvfz4sfai1BRKQpSREIABeP6sOw3K7MXlhKvWYJIiJfkDSBkJpizJxUwCfb9vHnj7eHXY6ISMxJmkAAuOz0fgzO6cKshWtw1yxBRCRSUgVCWmoKd0wo4MPNe3l1VXnY5YiIxJSkCgSAr53dn/49OnOfZgkiIkdJukBIT03hOxOG8e6nn/H22l1hlyMiEjOSLhAArho9gLxumdz3ypqwSxERiRlJGQid0lP51gXDKF5fwdL1FWGXIyISE5IyEACuGTuIXlkZzFqoWYKICCRxIHTOSOW284fyxpqdvPvp7rDLEREJXdIGAsD14wdzUpd0Zi8sDbsUEZHQJXUgdM1MY9p5+bzyyQ4+3Lwn7HJEREKV1IEAcOOXhpDdKU2zBBFJekkfCN06pXPLl4bw0sptrNq2L+xyRERCk/SBAHDrefl0zUhl9iLNEkQkeSkQgB5dMrjhnCH8z/tbWFteGXY5IiKhUCAEpp+fT2ZaCg8sWht2KSIioVAgBHplZXLduMH8fsVmPt11IOxyREQ6nAIhwowLhpKaYvzqNa0liEjyUSBEyOvWiauLBjJ/eRmbPzsYdjkiIh1KgdDItycMA+DB17SWICLJRYHQSP8enbny7AE8u2wTO/ZWhV2OiEiHUSA04fYJBdTVOw+9vi7sUkREOowCoQmDcrow5Yx+PFX8KbsqD4VdjohIh4gqEMzsYjNbZWalZvb9JrZnmtlzwfZiMxsStF9oZsvN7IPg46SIfUYH7aVmdp+ZWVsNqi3cPrGAqto65r65PuxSREQ6RIuBYGapwP3AJcBI4BozG9mo2zRgt7sXAL8A7g3adwKXuftpwE3AExH7/AqYARQGt4tPYBxtrqB3Fn99Wl8ef3sDnx2oDrscEZF2F80MYSxQ6u7r3L0aeBaY0qjPFGBecH8+MNnMzN3fdfctQftKoFMwm+gLdHP3xe7uwOPAFSc8mjY2c1IB+6vreOStDWGXIiLS7qIJhP7ApojHZUFbk33cvRbYA+Q06nMl8K67Hwr6l7VwTADMbIaZlZhZSXl5eRTltp2T+3TjolPzeOyt9eytqunQzy0i0tGiCYSmzu17a/qY2ak0nEb6ViuO2dDo/pC7F7l7UW5ubhTltq07JxWyt6qWJxZv7PDPLSLSkaIJhDJgYMTjAcCW5vqYWRrQHagIHg8Afgfc6O5rI/oPaOGYMWFU/+5MHJHL3DfWsf9QbdjliIi0m2gCYRlQaGb5ZpYBTAUWNOqzgIZFY4CrgIXu7mbWA3gB+IG7v3W4s7tvBfaZ2fjg6qIbgT+c4FjazZ2TC9l9oIanijVLEJHE1WIgBGsCM4GXgY+B5919pZndY2aXB90eBnLMrBS4Gzh8aepMoAD4oZmtCG69g23fAeYCpcBa4MW2GlRbO3vQSZxX0IuHXl9PVU1d2OWIiLQLa7jIJz4UFRV5SUlJKJ+7eN0urn5oCT++bCS3nJsfSg0iIsfDzJa7e1FL/fRK5SiNG5rD2PyePPjaOg7VapYgIolHgdAKd00qZNveKuYvL2u5s4hInFEgtMK5BTmcObAHv3p1LTV19WGXIyLSphQIrWBm3DW5gLLdB/ndu5vDLkdEpE0pEFpp4ojejOrfjQcWlVKrWYKIJBAFQiuZGTMnFrJh1wFe+GBr2OWIiLQZBcJx+MrIPEbkZTN7YSn19fFz2a6IyLEoEI5DSopxx6QC1uyo5KWV28IuR0SkTSgQjtNfn9aXobldmbWwlHh6cZ+ISHMUCMcpNcW4Y0IBH2/dy18+3hF2OSIiJ0yBcAKmnNmPQT27MGvhGs0SRCTuKRBOQFpqCrdPGMb7ZXt4fc3OsMsRETkhCoQT9PWzB9CveydmvaJZgojENwXCCcpIS+HbE4ZRsnE3i9ftCrscEZHjpkBoA98sGkjv7ExmvVIadikiIsdNgdAGOqWnMuOCoSxet4uSDRVhlyMiclwUCG3kunGDyemawX0LNUsQkfikQGgjnTNSmX7+UF5fXc57mz4LuxwRkVZTILShG84ZTPfO6czSLEFE4pACoQ1lZaZx67n5/OXj7Xy0ZW/Y5YiItIoCoY3dfO4QsjPTmL1oTdiliIi0igKhjXXvnM5NXxrCix9uY832fWGXIyISNQVCO7j1vHw6p6dy/yKtJYhI/FAgtIOeXTO4YfxgFry3hfU794ddjohIVBQI7WT6+UNJT03hAc0SRCROKBDaSW52JteMHcTv3t3MpooDYZcjItIiBUI7+vaXh5Fixq9eWxt2KSIiLVIgtKM+3TvxjaIBzC8pY+ueg2GXIyJyTAqEdvadCcOod+fB19aFXYqIyDEpENrZgJO68PWz+/PM0k/Zsa8q7HJERJqlQOgAt08ooKaunrlvrA+7FBGRZikQOsCQXl25/Ix+PLlkIxX7q8MuR0SkSQqEDjJzUgEHa+p4+E2tJYhIbIoqEMzsYjNbZWalZvb9JrZnmtlzwfZiMxsStOeY2SIzqzSz2Y32eTU45org1rstBhSrCnpnc+movsx7eyN7DtSEXY6IyBe0GAhmlgrcD1wCjASuMbORjbpNA3a7ewHwC+DeoL0K+CHwd80c/jp3PzO47TieAcSTmZMKqDxUy2Nvbwi7FBGRL4hmhjAWKHX3de5eDTwLTGnUZwowL7g/H5hsZubu+939TRqCIemd0rcbF47M45G31rOvSrMEEYkt0QRCf2BTxOOyoK3JPu5eC+wBcqI49qPB6aIfmpk11cHMZphZiZmVlJeXR3HI2HbnpAL2HKzhiSUbwy5FROQo0QRCU7+o/Tj6NHadu58GnB/cbmiqk7s/5O5F7l6Um5vbYrGx7vQBPfjy8FzmvrGeA9W1YZcjInJENIFQBgyMeDwA2NJcHzNLA7oDFcc6qLtvDj7uA56m4dRUUrhrcgEV+6t5uvjTsEsRETkimkBYBhSaWb6ZZQBTgQWN+iwAbgruXwUsdPdmZwhmlmZmvYL76cBXgQ9bW3y8Gj24J18alsODr6+jqqYu7HJERIAoAiFYE5gJvAx8DDzv7ivN7B4zuzzo9jCQY2alwN3AkUtTzWwD8HPgZjMrC65QygReNrP3gRXAZmBO2w0r9t05qZDyfYd4vmRTy51FRDqAHeMP+ZhTVFTkJSUlYZfRJtydbz64mM27D/Lq308kI02vERSR9mFmy929qKV++i0UEjNj5qRCtuyp4rfvlIVdjoiIAiFMFxT24owB3Xng1VJq6urDLkdEkpwCIURmxp2TCtlUcZA/rGh84ZaISMdSIIRs8im9Gdm3Gw8sKqWuPn7Wc0Qk8SgQQtYwSyhg3c79vPDB1rDLEZEkpkCIARed2ofC3lnMXriGes0SRCQkCoQYkJJizJxUwOrtlfzpo21hlyMiSUqBECO+eno/8nt1ZdbCUuLptSEikjgUCDEiNcW4fcIwVm7Zy6JVCf/WECISgxQIMeSKs/oz4KTO3PeKZgki0vEUCDEkPTWF2ycUsGLTZ7xZujPsckQkySgQYsyVo/vTt3snZr1SGnYpIpJkFAgxJjMtlW9dMJSlGypYsm5X2OWISBJRIMSgqWMH0Ssrk1kL14RdiogkEQVCDOqU3jBLeKt0F8s37g67HBFJEgqEGHXd+EH07JrBbM0SRKSDKBBiVJeMNKadl8+iVeV8ULYn7HJEJAkoEGLYjecMplunNK0liEiHUCDEsOxO6dxybj5/+mg7H2/dG3Y5IpLgFAgx7tZz88nKTGP2Ir0uQUTalwIhxnXvks6N5wzmjx9spXRHZdjliEgCUyDEgWnn5dMpLZUHNEsQkXakQIgDOVmZXD9+EH94bwsbd+0PuxwRSVAKhDhx2/lDSU0xHli0NuxSRCRBKRDiRO9unbhmzEB++04ZZbsPhF2OiCQgBUIc+daXh2EGD762LuxSRCQBKRDiSL8enblq9ECeK9nE9r1VYZcjIglGgRBnbp8wjLp61yxBRNqcAiHODOzZhSvO7M/TSzeys/JQ2OWISAJRIMShOyYOo7q2njlvaJYgIm1HgRCHhuZm8dXT+/HE4o3s3l8ddjkikiAUCHFq5qQCDlTX8ehb68MuRUQShAIhTg3Py+aSUX149O0N7K2qCbscEUkAUQWCmV1sZqvMrNTMvt/E9kwzey7YXmxmQ4L2HDNbZGaVZja70T6jzeyDYJ/7zMzaYkDJZOakAvZV1TLvrQ1hlyIiCaDFQDCzVOB+4BJgJHCNmY1s1G0asNvdC4BfAPcG7VXAD4G/a+LQvwJmAIXB7eLjGUAyO7Vfdyaf3JuH31pP5aHasMsRkTgXzQxhLFDq7uvcvRp4FpjSqM8UYF5wfz4w2czM3fe7+5s0BMMRZtYX6Obui93dgceBK05kIMnqzsmFfHaghieXbAy7FBGJc9EEQn9gU8TjsqCtyT7uXgvsAXJaOGZZC8cEwMxmmFmJmZWUl5dHUW5yOXNgD84v7MXcN9ZxsLou7HJEJI5FEwhNndv34+hzXP3d/SF3L3L3otzc3GMcMnndNbmQnZXVPLP007BLEZE4Fk0glAEDIx4PALY018fM0oDuQEULxxzQwjElSmOG9GT80J48+Ppaqmo0SxCR4xNNICwDCs0s38wygKnAgkZ9FgA3BfevAhYGawNNcvetwD4zGx9cXXQj8IdWVy9H3DWpkO17D/Gb5WUtdxYRaUKLgRCsCcwEXgY+Bp5395Vmdo+ZXR50exjIMbNS4G7gyKWpZrYB+Dlws5mVRVyh9B1gLlAKrAVebJshJadzhuVw9qAe/PrVtVTX1oddjojEITvGH/Ixp6ioyEtKSsIuI2YtWrWDWx5dxr1XnsbVYwaFXY6IxAgzW+7uRS310yuVE8iE4bmcPqA7D7y6lto6zRJEpHUUCAnEzJg5sYCNuw7w3+9rjV5EWkeBkGAuHJnHyX2ymb2wlLr6+DkdKCLhSwu7AGlbZsadkwq54+l3uPnRpXTrnB52SR0mIzWF6efnc2q/7mGXIhKXFAgJ6OJRffjKyDzWlley5bODYZfTYcr3HWLhJzt45rbxjOzXLexyROKOrjKShLGp4gDffHAxh2rreXbGeIbnZYddkkhM0FVGknQG9uzCM7eNJy3FuHZOMaU7KsMuSSSuKBAkoQzp1ZWnbxsPwLVzlrBh5/6QKxKJHwoESTgFvbN4avo4auuda+csYVPFgbBLEokLCgRJSCP6ZPPktHHsr67jmjlL2JxEi+six0uBIAlrZL9uPDltHHsO1nDtnCVs21PV8k4iSUyBIAnttAHdefzWseyqrObauUvYsU+hINIcBYIkvLMGncSjt4xh254qrptTzK7KQ2GXJBKTFAiSFMYM6cnDN41h0+4DXDe3mN37q8MuSSTmKBAkaZwzLIc5Nxaxbud+bnikmD0Ha8IuSSSmKBAkqZxfmMuDN4xm9bZKbnxkKfuqFAoihykQJOlMHNGb+687m5Wb93Dzo8vYf6g27JJEYoICQZLShSPzmHXNWazY9Bm3PraMg9V1YZckEjoFgiStS07ry8+/eQbLNlRw2+MlVNUoFCS5KRAkqU05sz8/veoM3lq7k289sZxDtQoFSV4KBEl6V44ewL9/7TReW13OHU+9Q3Wt3o9akpMCQQSYOnYQ/3rFKP7y8Q7ueuZdauoUCpJ8FAgigRvGD+ZHXx3JSyu3cffz7+k9qSXp6C00RSLcel4+NXX1/PuLn5CeYvz0G2eQmmJhlyXSIRQIIo1868vDqK6t52d/Xk16agr//vXTSFEoSBJQIIg04c7JhdTU1XPfwlLS04x/nTIKM4WCJDYFgkgz/ubC4VTXOb9+bS3pqSn86KsjFQqS0BQIIs0wM7538Qiqa+t55K31ZKSm8P1LTlYoSMJSIIgcg5nxw6+eQk1dPQ++vo6MtBT+9isjwi5LpF0oEERaYGb8y+WnUltfz6yFpaSnpnDX5MKwyxJpcwoEkSikpBj/dsVpVNc6Pw+uPvrOhGFhlyXSphQIIlFKSTH+86rTqa2v596XPiE91Zh+/tCwyxJpM1G9UtnMLjazVWZWambfb2J7ppk9F2wvNrMhEdt+ELSvMrOLIto3mNkHZrbCzEraYjAi7S01xfjZN87g0tP68JMXPubxxRvCLkmkzbQ4QzCzVOB+4EKgDFhmZgvc/aOIbtOA3e5eYGZTgXuBq81sJDAVOBXoB/zFzIa7++F/KTnR3Xe24XhE2l1aagr/NfUsaure4Ud/WElaSgrXjhsUdlkiJyyaGcJYoNTd17l7NfAsMKVRnynAvOD+fGCyNVybNwV41t0Puft6oDQ4nkhcS09NYfa1ZzFxRC7/+PsP+E3JprBLEjlh0QRCfyDy2V4WtDXZx91rgT1ATgv7OvAnM1tuZjOa++RmNsPMSsyspLy8PIpyRTpGZloqv7p+NOcV9OIffvs+v393c9gliZyQaAKhqVfhNP43kM31Oda+57r72cAlwB1mdkFTn9zdH3L3Incvys3NjaJckY7TKT2Vh24oYnx+Dnc/v4IX3t8adkkixy2aQCgDBkY8HgBsaa6PmaUB3YGKY+3r7oc/7gB+h04lSZzqnJHK3JuKGD34JP7Xs+/y8sptYZckclyiCYRlQKGZ5ZtZBg2LxAsa9VkA3BTcvwpY6O4etE8NrkLKBwqBpWbW1cyyAcysK/AV4MMTH45IOLpmpvHoLWM5bUB3Zj79Dgs/2R52SSKt1mIgBGsCM4GXgY+B5919pZndY2aXB90eBnLMrBS4G/h+sO9K4HngI+Al4I7gCqM84E0zew9YCrzg7i+17dBEOlZWZhqP3TKWk/t049tPvMPrq7XmJfHFGv6Qjw9FRUVeUqKXLEhs++xANdfMKWZdeSWP3jyGLxX0CrskSXJmttzdi1rqp7fQFGljPbpk8NT0cQzO6cK0eSUsXV8RdkkiUVEgiLSDnl0zeGr6ePr16MQtjy5l+cbdYZck0iIFgkg7yc3O5OnbxpObncnNjyzl/bLPwi5J5JgUCCLtKK9bJ56+bTw9uqZz/dxiPty8J+ySRJqlQBBpZ/16dObp6ePJykzjhoeL+WTb3rBLEmmSAkGkAwzs2YVnZownIy2F6+cWU7pjX9gliXyBAkGkgwzO6crTt40H7MhlqSKxRIEg0oGG5WbxzG3jqK93rp1TzKe7DoRdksgRCgSRDlaYl82T08dRVVvHNXOWULZboSCxQYEgEoJT+nbjyWnj2FdVw7Vzitm652DYJYkoEETCMqp/dx6fNo6K/dVcN6eYHXurwi5JkpwCQSREZw7swbxbx7BtbxXXzi1mZ+WhsEuSJKZAEAnZ6ME9efTmMZTtPsD1c4vZvb867JIkSSkQRGLAuKE5PHzTGNbv3M/1Dxez50BN2CVJElIgiMSIcwt68eANo1mzvZIbHylmb5VCQTqWAkEkhkwY0ZsHrjublVv2cvMjS6k8VBt2SZJEFAgiMeavRuYx+9qzeK9sD7c+towD1QoF6RgKBJEYdPGovvzy6jMp2VDB9HklVNXUhV2SJAEFgkiMuuyMfvzsm2eweN0uZjyxXKEg7U6BIBLDvnbWAO79+um8vrqcO556h+ra+rBLkgSmQBCJcd8cM5CfXDGKVz7ZwZ3PvENNnUJB2ocCQSQOXD9+MD++bCQvr9zOd59bQa1CQdpBWtgFiEh0bjk3n9o659/++DEZqSn832+cQWqKhV2WJBAFgkgcue2CoVTX1fPTl1eRlmLce+XppCgUpI0oEETizB0TC6iuree/XllDWqrxt18ZQa+szLDLkgSgQBCJQ9/9q0Jq6up54NW1PLN0Ez27ZjA8L4sRedkM75PN8LxshvfOpnuX9LBLlTiiQBCJQ2bG3180ggkjevPh5j2s3r6PVdv3MX95GfurP3+9Qp9unRjeJ5sReVkU5mUzIi+bwrwsumToR1++SM8KkThlZozN78nY/J5H2tydLXuqWL2tISAOf3x83S4ORbyGYVDPLgzPy2J4XjYjghnF0NyuZKalhjEUiREKBJEEYmb079GZ/j06M/Hk3kfa6+qdTysOsGrbPlZv//z26qpyausdgNQUY0hOlyMBcfg2JKcLaam6Qj0ZKBBEkkBqipHfqyv5vbpy8ag+R9qra+tZv3P/kYBYtW0fH23Zy4sfbsMbcoKM1BSG9c466rTTiD7Z9O/RWVc4JRgFgkgSy0hLYUSfhl/wkQ5W17G2vPLIjGLV9n0s27Cb36/YcqRPl4xUCnsffdppRJ9semdnYqagiEcKBBH5gs4ZqYzq351R/bsf1b63qoY12yuPOu20aFU5v1ledqRPt05pRwVEYe+Gjz27ZnT0MKSVFAgiErVundIZPfgkRg8+6aj2iv3VR512Wr19H//93haeKv78vRx6ZWUyok/WkYBoWKPIIruTLo2NFVEFgpldDPwXkArMdff/aLQ9E3gcGA3sAq529w3Bth8A04A64C53fzmaY4pI/OjZNYPxQ3MYPzTnSJu7s2Pfoc9PO23bx+odlTxfsokDEZfG9ut++NLYzxeyC3pn0TlDVzx1tBYDwcxSgfuBC4EyYJmZLXD3jyK6TQN2u3uBmU0F7gWuNrORwFTgVKAf8BczGx7s09IxRSSOmRl53TqR160TFwzPPdJeX+9s/uxgEBCHL42t5O3SXVQH/7TPDAb37PL5aadgMTu/V1cy0nTFU3uJZoYwFih193UAZvYsMAWI/OU9Bfjn4P58YLY1rCpNAZ5190PAejMrDY5HFMcUkQSUkmIM7NmFgT278Fcj846019bVs7HiwOevodi+j9XbK3nlkx3UBZfGpgX7piXh1U3/c9d57f46kWgCoT+wKeJxGTCuuT7uXmtme4CcoH1Jo337B/dbOiYAZjYDmAEwaNCgKMoVkXiUlprCsNwshuVmcclpfY+0H6qtY13555fGbth1AD98TWwSMdo/BKMJhKaqaPzdaK5Pc+1Nzfma/A67+0PAQwBFRUXJ9ywQSXKZaamc0rcbp/TtFnYpCS+ak3FlwMCIxwOALc31MbM0oDtQcYx9ozmmiIh0oGgCYRlQaGb5ZpZBwyLxgkZ9FgA3BfevAhZ6w5xuATDVzDLNLB8oBJZGeUwREelALZ4yCtYEZgIv03CJ6CPuvtLM7gFK3H0B8DDwRLBoXEHDL3iCfs/TsFhcC9zh7nUATR2z7YcnIiLRsnhanCkqKvKSkpKwyxARiStmttzdi1rqpwt6RUQEUCCIiEhAgSAiIoACQUREAnG1qGxm5cDG49y9F7CzDcuJBxpzcki2MSfbeOHExzw7N5PeAAAEU0lEQVTY3XNb6hRXgXAizKwkmlX2RKIxJ4dkG3OyjRc6bsw6ZSQiIoACQUREAskUCA+FXUAINObkkGxjTrbxQgeNOWnWEERE5NiSaYYgIiLHoEAQEREgQQPBzB4zs/VmtiK4nRm0m5ndZ2alZva+mZ0ddq0nwsxmBmNxM+sV0d7sOM3sJjNbE9xuavrIscvMnjKzVWb2oZk9YmbpQXvCjvkwM5tlZpURjzPN7LlgzMVmNiRi2w+C9lVmdlEY9Z4IM5tsZu8EP79vmllB0J4wY27u5zdi+xgzqzOzqyLamnwum9loM/sgON59wVsYt567J9wNeAy4qon2S4EXaXgnt/FAcdi1nuA4zwKGABuAXi2NE+gJrAs+nhTcPynscbRyzJcG4zLgGeA7iT7mYBxFwBNAZUTb7cCvg/tTgeeC+yOB94BMIB9YC6SGPYZWjnc1cErEOB9LtDE39/MbbEsFFgJ/PPy77FjPZRreZ+ac4Pn/InDJ8dSUkDOEY5gCPO4NlgA9zKxvSzvFKnd/1903NLGpuXFeBPzZ3SvcfTfwZ+Dijqv4xLn7H4NxOQ0/BAOCTQk7ZjNLBX4K/EOjTVOAecH9+cDk4C/DKcCz7n7I3dcDpcDYjqq3jThw+D0zu/P5OyomzJiP8fMLcCfwW2BHRFuTz+Xged7N3RcHPxePA1ccT02JHAj/Fpw6+IWZZQZt/YFNEX3KgrZE09w4E2b8wamiG4CXgqZEHvNMYIG7b23UfmRs7l4L7AFySIwxTwf+aGZlNHyf/yNoT+QxA2Bm/YGvAb9utOlYz/GyJtpbLVED4QfAycAYGqZX3wvamzqvlojX3TY3zkQa/wPA6+7+RvA4IcdsZv2AbwCzmtrcRFvcjznwN8Cl7j4AeBT4edCeyGM+7JfA9zx4d8kI7T72hAwEd98anDo4RMOT6fDUsQwYGNF1AJ9PRRNJc+NMiPGb2Y+BXODuiOZEHfNZQAFQamYbgC7BW9VCxNjMLI2GUysVxPmYzSwXOMPdi4Om54AvBfcTcsyNFAHPBt/vq4AHzOwKjv0cH9BEe+uFvbDSTos1fYOPRkPa/kfw+K85euFxadi1ttF4N3D0onKT46RhtrSehgWpk4L7PcOuv5VjnQ68DXRu1J6wY240zshF5Ts4eoH1+eD+qRy9wLqOGF9gbTTGNBr+s+fw4PE04LeJOubGP7+Ntj3G0YvKTT6XgWXB8/7wovKlx1VL2F+MNv7C/hHoR8Pq/AfAh8CTQFaw3YD7abgC4QOgKOyaT3Ccd9Hw10EtDX8RzG1pnMCtNCy4lQK3hD2W4xhzbTCuFcHtR4k+5kZtkYHQCfhNMK6lwNCIbf8YfC1WcZxXnIT8ff5a8H18D3j18NgSYcwt/fw26nskEILHTT6XaZhVfBiMfzbBf6Fo7U3/ukJERIAEXUMQEZHWUyCIiAigQBARkYACQUREAAWCiIgEFAgiIgIoEEREJPD/ATWFaq4dqHIYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "tau = 8.\n",
    "\n",
    "def readMatrix(file):\n",
    "    fd = open(file, 'r')\n",
    "    hdr = fd.readline()\n",
    "    rows, cols = [int(s) for s in fd.readline().strip().split()]\n",
    "    tokens = fd.readline().strip().split()\n",
    "    matrix = np.zeros((rows, cols))\n",
    "    Y = []\n",
    "    for i, line in enumerate(fd):\n",
    "        nums = [int(x) for x in line.strip().split()]\n",
    "        Y.append(nums[0])\n",
    "        kv = np.array(nums[1:])\n",
    "        k = np.cumsum(kv[:-1:2])\n",
    "        v = kv[1::2]\n",
    "        matrix[i, k] = v\n",
    "    #化为正负1\n",
    "    category = (np.array(Y) * 2) - 1\n",
    "    return matrix, tokens, category\n",
    "\n",
    "def svm_train(matrix, category):\n",
    "    state = {}\n",
    "    M, N = matrix.shape\n",
    "    #####################\n",
    "    #大于0的化为1\n",
    "    matrix = 1.0 * (matrix > 0)\n",
    "    '''\n",
    "    #构造kernel矩阵\n",
    "    d1 = np.sum(matrix ** 2, axis=1).reshape(-1, 1)\n",
    "    d2 = d1.T\n",
    "    squared = matrix.dot(matrix.T)\n",
    "    dist = d1 + d2 - 2 * squared\n",
    "    k = np.exp(- dist / (2 * (tau ** 2)))\n",
    "    '''\n",
    "    gram = matrix.dot(matrix.T)\n",
    "    squared = np.sum(matrix*matrix, axis=1)\n",
    "    k = np.exp(-(squared.reshape((-1,1)) + squared.reshape((1,-1)) - 2 * gram) / (2 * (tau ** 2)))\n",
    "    \n",
    "    #初始化向量\n",
    "    alpha = np.zeros(M)\n",
    "    #循环次数\n",
    "    n = 40\n",
    "    #系数\n",
    "    L = 1. / (64 * M)\n",
    "    #平均值\n",
    "    alpha_avg = np.zeros(M)\n",
    "    \n",
    "    for j in range(n * M):\n",
    "        #随机取一个更新的维度\n",
    "        i = int(np.random.rand() * M)\n",
    "        #计算函数间隔\n",
    "        margin = category[i] * (k[i, :].dot(alpha))\n",
    "        #只改变第i个分量\n",
    "        grad = M * L * k[:, i] * alpha[i]\n",
    "        if(margin < 1):\n",
    "            grad -= category[i] * k[:, i]\n",
    "        alpha -= grad / ((np.sqrt(j+1)))\n",
    "        alpha_avg += alpha\n",
    "    \n",
    "    alpha_avg /= (n * M)\n",
    "    \n",
    "    state['alpha'] = alpha\n",
    "    state['alpha_avg'] = alpha_avg\n",
    "    state['Xtrain'] = matrix\n",
    "    state['Sqtrain'] = squared\n",
    "    \n",
    "    ####################\n",
    "    return state\n",
    "\n",
    "\n",
    "def svm_test(matrix, state):\n",
    "    M, N = matrix.shape\n",
    "    output = np.zeros(M)\n",
    "    ###################\n",
    "    Xtrain = state['Xtrain']\n",
    "    Sqtrain = state['Sqtrain']\n",
    "    #大于0的化为1\n",
    "    matrix = 1.0 * (matrix > 0)\n",
    "    #做测试集的kernel\n",
    "    gram = matrix.dot(Xtrain.T)\n",
    "    squared = np.sum(matrix * matrix, axis=1)\n",
    "    k = np.exp(-(squared.reshape((-1, 1)) + Sqtrain.reshape((1, -1)) - 2 * gram) / (2 * (tau**2)))\n",
    "    #读取alpha\n",
    "    alpha_avg = state['alpha_avg']\n",
    "    #预测\n",
    "    pred = k.dot(alpha_avg)\n",
    "    output = np.sign(pred)\n",
    "    \n",
    "    ###################\n",
    "    return output\n",
    "\n",
    "def evaluate(output, label):\n",
    "    error = (output != label).sum() * 1. / len(output)\n",
    "    print('Error: %1.4f' % error)\n",
    "    return error\n",
    "\n",
    "def svm(file):\n",
    "    trainMatrix, tokenlist, trainCategory = readMatrix(file)\n",
    "    testMatrix, tokenlist, testCategory = readMatrix('MATRIX.TEST')\n",
    "    \n",
    "    state = svm_train(trainMatrix, trainCategory)\n",
    "    output = svm_test(testMatrix, state)\n",
    "    \n",
    "    return evaluate(output, testCategory)\n",
    "\n",
    "\n",
    "trainMatrix, tokenlist, trainCategory = readMatrix('MATRIX.TRAIN.400')\n",
    "testMatrix, tokenlist, testCategory = readMatrix('MATRIX.TEST')\n",
    "\n",
    "state = svm_train(trainMatrix, trainCategory)\n",
    "output = svm_test(testMatrix, state)\n",
    "\n",
    "evaluate(output, testCategory)\n",
    "\n",
    "size = ['.50','.100','.200','.400','.800','.1400']\n",
    "size1 = [50, 100, 200, 400, 800, 1400]\n",
    "train = \"MATRIX.TRAIN\"\n",
    "error = []\n",
    "for i in size:\n",
    "    file = train+i\n",
    "    error.append(svm(file))\n",
    "    \n",
    "plt.plot(size, error)\n",
    "plt.title(\"error VS szie\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error: 0.0025\n",
      "Error: 0.0163\n",
      "Error: 0.0075\n",
      "Error: 0.0025\n",
      "Error: 0.0025\n",
      "Error: 0.0000\n",
      "Error: 0.0000\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'error VS szie')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEICAYAAABfz4NwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XucVXW9//HXmxm5aAmIk8oMOiAYoomXiV/mJYvTEUwZLc4vPKfS0qiUrOMpxTpZkaeyi5eOduFoZeURiC5Sophplh1DBu/cdASNEdQxTNSOEPI5f+zv0Ha3Z2bPdc/e+/18PObh2t/1XWt/vzgz773W+sxaigjMzMwGFXsAZmY2MDgQzMwMcCCYmVniQDAzM8CBYGZmiQPBzMwAB4JZSZK0v6QXJVUVeyxWPhwIZl0k6WhJL0l6bZ5190mak5bPkrRW0guSnpZ0U75tuiMi/hgRr4mIV3pjf2bgQLAyJ6m6kLau7CMi7gZagHfl9DsUmATcIOktwBeB0yPitcDBwKKujd6sfzkQrORIGi3pJ5JaJW2QdF7Wus9JWizpR5K2Ame20zZE0hWSNqWvKyQNSfs4QVKLpAslPQV8L88wrgPel9P2PuCmiPgT8Ebg7oi4DyAitkTEdRHxQjtzOlPS+nQ0sUHSv6T2B9KpobavSOOrT8vVqd9wSddK2izpSUmX+HSSdZUDwUqKpEHAL4AHgFpgKvBxSSdmdWsEFgMjgOvbafs08CbgcGAyMAX496x97AvsBRwAzM4zlB8Cx0naP2tc/wz8IK1fDpwo6fOSjmkLm3bmtAfwDWB6Opp4M3A/QERMTqeGXgOcD6wD7s2zm+uAHcB44AjgH4Gz23tPs3wcCFZq3gjURMS8iNgeEeuB/wJmZfW5OyJ+HhE7I+J/22n7F2BeRDwTEa3A54H3Zu1jJ/DZiNiWtY9dImIjcCfwntQ0FRgK3JTW/w54J3BkavuTpMs6+NS+EzhU0rCI2BwRq7JXSjoWuASYERFbc9btA0wHPh4RL0XEM8DlOf8mZp1yIFipOQAYLenPbV/Ap4B9svpszLNdbtto4Ims10+ktjatEfFyJ2PJPm30XuC/I+KvbSsj4uaIOIXMkUYjcCZ5PrVHxEvAu4EPA5vTxeeJbesljSFz/eGMiHgkzzgOAHZL27b9m3wHeF0n4zd7FQeClZqNwIaIGJH19dqIOCmrT75b+Oa2bSLzi7TN/qmto33k+ilQK+mtZI4GfpCvUzoq+TVwO3BoO32WRcTbgf2AtWSOepA0DPg5cEVE3NzOODYC24C9s/5N9oyIQwqYg9kuDgQrNfcAW9MF32GSqiQdKumNXdzPDcC/S6qRtDdwMfCjruwgfbJfTOai8xMR0dS2TlKjpFmSRipjCvAW4A+5+5G0j6QZ6VrCNuBFoK2c9LvA2oj4Sgfj2AzcCnxd0p6SBkk6MFU6mRXMgWAlJdXdn0LmYvAG4FngGmB4F3d1CdAEPAg8ROZC7SXdGNJ1ZI40co8OngM+CDwKbCUTNl+NiOv5e4OAfyNzhLKFTHCck9bNAk7LqTQ6Ls8+3gcMBlan915M5mjDrGDyA3LMzAx8hGBmZokDwczMAAeCmZklDgQzMwOgSzf5Kra999476uvriz0MM7OSsnLlymcjoqazfiUVCPX19TQ1NXXe0czMdpH0ROe9fMrIzMwSB4KZmQEOBDMzSxwIZmYGOBDMzCxxIJiZGeBAMDOzpOwDYefOYOGKP3LLw08VeyhmZgNaQYEgaZqkdZKaJc3Ns36IpIVp/XJJ9al9lKQ70j3cr8rZZrCk+ZIekbRW0rt6Y0L5/OgPf+RzS1bxl+07+uotzMxKXqeBkB4KfjWZh3hPAk6XNCmn21nAcxExnszDvS9N7S8DnwE+kWfXnwaeiYiD0n7v7NYMOjFokLj4lEk8tfVlvnPn+r54CzOzslDIEcIUoDki1kfEdmABmQeGZ2sk8+QoyDypaaokRcRLEXEXmWDI9QHgS7DrmbPPdmsGBXhj/V6847D9+M5vH2PTn/+3r97GzKykFRIItWQe4t2mJbXl7RMRO4DngVHt7VDSiLT4BUn3SvqxpH3a6TtbUpOkptbW1gKGm99F0yeyM+DSW9Z2ex9mZuWskEBQnrbc524W0idbNVAH/D4ijgTuBr6Wr2NEzI+IhohoqKnp9GZ97aobuTuzjxvHjfdvYuUTz3V7P2Zm5aqQQGgBxmS9riPzMPC8fSRVk3ng+ZYO9vkn4C/Az9LrHwNHFjCWHvnICQfyutcO4Qu/XM3OnX6WtJlZtkICYQUwQdJYSYOBWcCSnD5LgDPS8kzg9oho9zduWvcL4ITUNBVY3YVxd8seQ6q5YNpE7t/4Z2584Mm+fjszs5LSaSCkawJzgGXAGmBRRKySNE/SjNTtWmCUpGbgfGBXaaqkx4HLgDMltWRVKF0IfE7Sg8B7gX/rpTl16J1H1HJY3XAuvXmdy1DNzLKogw/yA05DQ0P0xgNymh7fwsxv3815Uydw/tsP6oWRmZkNXJJWRkRDZ/3K/i+V82mo34uTD9uP79z5GE+6DNXMDKjQQACYO30iAF9xGaqZGVDBgVA3cndmH+8yVDOzNhUbCAAffsuB7LPnEOa5DNXMrLIDYY8h1Vxw4kQe2Phnfn6/y1DNrLJVdCAAnHZELZPrhnPpLWtdhmpmFa3iA6HtbqhPb93Gt303VDOrYBUfCABHHbAXp0we7TJUM6toDoSkrQz10ptdhmpmlcmBkNSOGMaHjh/Hkgc2sfKJju7LZ2ZWnhwIWT7UVob6C5ehmlnlcSBk2WNINRdOm8gDLc+7DNXMKo4DIceph9cyecwILr1lLS9tcxmqmVUOB0KOQYPExSdnylC/c+djxR6OmVm/cSDkcdQBI5kxeTTf+e16Wp77S7GHY2bWLwoKBEnTJK2T1Cxpbp71QyQtTOuXS6pP7aMk3SHpRUlXtbPvJZIe7skk+sKF0yciwaW3rCv2UMzM+kWngSCpCrgamA5MAk7PeupZm7OA5yJiPHA5cGlqfxn4DPCJdvb9TuDF7g29b9WOGMbs4w/kFy5DNbMKUcgRwhSgOSLWR8R2YAHQmNOnEbguLS8GpkpSRLwUEXeRCYZXkfQaMo/bvKTbo+9jH37LOPbdcyifdxmqmVWAQgKhFtiY9bolteXtk57B/DwwqpP9fgH4OjBgT9LvPriaC6e/ngdbnudn97kM1czKWyGBoDxtuR+XC+nzt87S4cD4iPhZp28uzZbUJKmptbW1s+69rnGyy1DNrDIUEggtwJis13XApvb6SKoGhgMdnXg/GjhK0uPAXcBBkn6Tr2NEzI+IhohoqKmpKWC4vWvQIPHZUybxzAvb+LbLUM2sjBUSCCuACZLGShoMzAKW5PRZApyRlmcCt0dEu0cIEfGtiBgdEfXAscAjEXFCVwffX47cfySNh49mvstQzayMdRoI6ZrAHGAZsAZYFBGrJM2TNCN1uxYYJamZzIXiXaWp6SjgMuBMSS15KpRKwoXTMmWoX/bdUM2sTKmDD/IDTkNDQzQ1NRXt/S//1SNc+etHWfzho2mo36to4zAz6wpJKyOiobN+/kvlLviQy1DNrIw5ELpg98HVzJ0+kYeefJ6fugzVzMqMA6GLZkwezeFjRvAVl6GaWZlxIHTRoEHi4lSG+q3fuAzVzMqHA6Ebjtx/JKcePpr5v3MZqpmVDwdCN10wbSKDBF9yGaqZlQkHQjeNHjGMD7/lQG56cDMrHvfdUM2s9DkQeuBDxx/IfsOHMs9lqGZWBhwIPTBscBUXTsuUof7k3pZiD8fMrEccCD3UePhojth/BF9Zts5lqGZW0hwIPSSJi0+eROsL2/jmb5qLPRwzs25zIPSCI/YfyWlH1PJfv9vAxi0uQzWz0uRA6CUXTHs9g3w3VDMrYQ6EXrLf8FSG+tBm7tngMlQzKz0OhF60qwz1l6tchmpmJceB0IuGDa5i7vSJPPzkVha7DNXMSkxBgSBpmqR1kpolzc2zfoikhWn9ckn1qX2UpDskvSjpqqz+u0u6SdJaSaskfbm3JlRsMyZnylC/umwdL7oM1cxKSKeBIKkKuBqYDkwCTs/zGMyzgOciYjxwOXBpan8Z+AzwiTy7/lpETASOAI6RNL17UxhYJPHZUw7JlKHe4TJUMysdhRwhTAGaI2J9RGwHFgCNOX0agevS8mJgqiRFxEsRcReZYNglIv4SEXek5e3AvUBdD+YxoBw+ZgTvPKKWa+5yGaqZlY5CAqEW2Jj1uiW15e0TETuA54FRhQxA0gjgFODX7ayfLalJUlNra2shuxwQPjnt9VRJLkM1s5JRSCAoT1tuCU0hff5+x1I1cAPwjYhYn69PRMyPiIaIaKipqel0sANFdhnq8vV/KvZwzMw6VUggtABjsl7XAZva65N+yQ8HCinGnw88GhFXFNC35Mw+fhyjhw9l3i9X84rLUM1sgCskEFYAEySNlTQYmAUsyemzBDgjLc8Ebo+IDn8DSrqETHB8vGtDLh3DBldx4fSJrNq0lZ+sdBmqmQ1snQZCuiYwB1gGrAEWRcQqSfMkzUjdrgVGSWoGzgd2laZKehy4DDhTUoukSZLqgE+TqVq6V9L9ks7uzYkNFDMmj+bIdDdUl6Ga2UCmTj7IDygNDQ3R1NRU7GF02f0b/8ypV/+ec044kAumTSz2cMyswkhaGRENnfXzXyr3A5ehmlkpcCD0kwumTaRK4ks3ryn2UMzM8nIg9JN9hw/lIyccyNKHnuIPLkM1swHIgdCP2spQv+AyVDMbgBwI/WjoblXMPelgl6Ga2YDkQOhnpxy2H0cdMJKvLFvHCy//tdjDMTPbxYHQzyRx8cmTePbFbXzzN48VezhmZrs4EIpg8pgRvPPIWq793Qb++CeXoZrZwOBAKJILp02kapDLUM1s4HAgFMk+ew7lnBMO5OaHXYZqZgODA6GIPnj8OGpHDGPeL1yGambF50AooqG7VTF3+kRWb97K4pUbO9/AzKwPORCK7OTD9qPhgJF81WWoZlZkDoQik8TFp0zi2Re3c/UdLkM1s+JxIAwAh9WN4F1H1vHduzbwxJ9eKvZwzKxCFRQIkqZJWiepWdLcPOuHSFqY1i+XVJ/aR0m6Q9KLkq7K2eYoSQ+lbb4hKd9zmSvGBdNeT3WV+NLStcUeiplVqE4DQVIVcDUwncwTzk6XNCmn21nAcxExHrgcuDS1vwx8BvhEnl1/C5gNTEhf07ozgXLRVoZ6y6qnuPsxl6GaWf8r5AhhCtAcEesjYjuwAGjM6dMIXJeWFwNTJSkiXoqIu8gEwy6S9gP2jIi707OXfwCc2pOJlIOzj0tlqL4bqpkVQSGBUAtk10S2pLa8fdIzmJ8HRnWyz+zbfebbJwCSZktqktTU2tpawHBL19DdqrjopIms2byVHze5DNXM+lchgZDv3H7ux9dC+nSrf0TMj4iGiGioqanpYJfl4R1v2I831o/ka7e6DNXM+lchgdACjMl6XQdsaq+PpGpgOLClk33WdbLPipS5G+oh/Oml7Vx1R3Oxh2NmFaSQQFgBTJA0VtJgYBawJKfPEuCMtDwTuD1dG8grIjYDL0h6U6oueh9wY5dHX6beUDecdx1Zx/fuetxlqGbWbzoNhHRNYA6wDFgDLIqIVZLmSZqRul0LjJLUDJwP7CpNlfQ4cBlwpqSWrAqljwDXAM3AY8DNvTOl8vDJEzNlqF9c6ruhmln/qC6kU0QsBZbmtF2ctfwy8E/tbFvfTnsTcGihA600++w5lHPfOp6vLlvH/zz2LG8+cO9iD8nMypz/UnkAO+vYsdSOGMYXfrnGZahm1uccCAPY0N2q+NRJB7Nm81YWuQzVzPqYA2GAO+kN+2bKUJetY6vLUM2sDzkQBri2MtQtf9nO1be7DNXM+o4DoQS8oW44M4+s47u/38Djz7oM1cz6hgOhRHzyxNczuGqQy1DNrM84EErE6/YcyjlvHc+tq5/mf5qfLfZwzKwMORBKyFnHjqVupO+GamZ9w4FQQtrKUNc+9QILV7gM1cx6lwOhxEw/dF+m1O/F1291GaqZ9S4HQomRxMWnTGLLX7ZzlctQzawXORBK0KG1w/mno+r4nstQzawXORBK1CdchmpmvcyBUKJe99qhnPs2l6GaWe9xIJSwDxzjMlQz6z0FBYKkaZLWSWqWNDfP+iGSFqb1yyXVZ627KLWvk3RiVvu/Slol6WFJN0ga2hsTqiTZZagLVvyx2MMxsxLXaSBIqgKuBqYDk4DTs5561uYs4LmIGA9cDlyatp1E5pGbhwDTgG9KqpJUC5wHNETEoUBV6mddNP3QfZkydi++fusjLkM1sx4p5AhhCtAcEesjYjuwAGjM6dMIXJeWFwNT07OSG4EFEbEtIjaQeVzmlNSvGhgmqRrYHdjUs6lUpszdUCfxnMtQzayHCgmEWiD7z2JbUlvePukZzM8Do9rbNiKeBL4G/BHYDDwfEbd2ZwKWKUP9/0eN4Xu/38AGl6GaWTcVEgjK05Z7BbO9PnnbJY0kc/QwFhgN7CHpPXnfXJotqUlSU2trawHDrUz/duJBLkM1sx4pJBBagDFZr+v4+9M7u/qkU0DDgS0dbPsPwIaIaI2IvwI/Bd6c780jYn5ENEREQ01NTQHDrUxtZai/Wv00v3cZqpl1QyGBsAKYIGmspMFkLv4uyemzBDgjLc8Ebo+ISO2zUhXSWGACcA+ZU0VvkrR7utYwFfBH2x76wDFjGbPXMOb9YjU7XtlZ7OGYWYnpNBDSNYE5wDIyv7QXRcQqSfMkzUjdrgVGSWoGzgfmpm1XAYuA1cAtwLkR8UpELCdz8fle4KE0jvm9OrMKNHS3Kj41/WDWPf0CC5t8N1Qz6xplPsiXhoaGhmhqair2MAa0iGDW/D/w6DMvcscnTmD4sN2KPSQzKzJJKyOiobN+/kvlMiOJz6Qy1It++qBPHZlZwRwIZejQ2uF8avrBLH3oKT6+8H6HgpkVpLrYA7C+8cHjx7Ezgi/dvBaAK959ONVVzn8za58DoYx96C0HAjgUzKwgDoQy51Aws0I5ECrAh95yIBJ8celaArjSoWBmeTgQKsTs4zNHCl9cmjlScCiYWS4HQgVxKJhZRxwIFcahYGbtcSBUIIeCmeXjQKhQs48/ECH+Y+kaCLhylkPBrNI5ECrYB48fB5AJBRwKZpXOgVDhHApm1saBYA4FMwMcCJbkhsIVsw5nN4eCWUVxINguHzx+HBJccpNDwawSFfTTLmmapHWSmiXNzbN+iKSFaf1ySfVZ6y5K7esknZjVPkLSYklrJa2RdHRvTMh65uzjxvHv7ziYmx7azMcX3M9ffetss4rR6RGCpCrgauDtQAuwQtKSiFid1e0s4LmIGC9pFnAp8G5Jk8g8g/kQYDRwm6SDIuIV4ErgloiYmZ7VvHuvzsy67ezjMqePfKRgVlkKOWU0BWiOiPUAkhYAjWSek9ymEfhcWl4MXCVJqX1BRGwDNqRnLk+RtAo4HjgTICK2A9t7PBvrNdmhEARXzjrCoWBW5gr5Ca8Fsp/Y3pLa8vaJiB3A88CoDrYdB7QC35N0n6RrJO2R780lzZbUJKmptbW1gOFab2k7fbT0oaf42IL7fPrIrMwVEgjK0xYF9mmvvRo4EvhWRBwBvAT83bUJgIiYHxENEdFQU1NTwHCtNzkUzCpHIYHQAozJel0HbGqvj6RqYDiwpYNtW4CWiFie2heTCQgbgBwKZpWhkEBYAUyQNDZd/J0FLMnpswQ4Iy3PBG6PiEjts1IV0lhgAnBPRDwFbJT0+rTNVF59TcIGGIeCWfnr9KJyROyQNAdYBlQB342IVZLmAU0RsQS4Fvhhumi8hUxokPotIvPLfgdwbqowAvgocH0KmfXA+3t5btbLXnWhOe7jG6f7QrNZOVHmg3xpaGhoiKampmIPo+Jd87v1XHLTGqYfuq9DwawESFoZEQ2d9fNPsnXZ2ceN4zMnT+Lmh5/ivBt8+sisXDgQrFvOOnasQ8GszDgQrNscCmblxYFgPeJQMCsfDgTrsexQ+Oh/OxTMSpUDwXrFWceO5eKTJ3HLKoeCWalyIFiv+YBDwaykORCsVzkUzEqXA8F6XXYozPnvex0KZiXCgWB94gPHjuWzp0xi2aqnHQpmJcKBYH3m/cc4FMxKiQPB+pRDwax0OBCsz+WGwvYdDgWzgciBYP0iOxQ+eoNDwWwgciBYv3n/MWP5nEPBbMAqKBAkTZO0TlKzpL979nF6ItrCtH65pPqsdRel9nWSTszZrkrSfZJ+2dOJWGk406FgNmB1GgiSqoCrgenAJOB0SZNyup0FPBcR44HLgUvTtpPIPD3tEGAa8M20vzYfA9b0dBJWWhwKZgNTIUcIU4DmiFgfEduBBUBjTp9G4Lq0vBiYKkmpfUFEbIuIDUBz2h+S6oB3ANf0fBpWarJDwReazQaGQgKhFtiY9bolteXtExE7gOeBUZ1sewVwAdDhbwJJsyU1SWpqbW0tYLhWKtpC4dbVDgWzgaCQQFCettwHMbfXJ2+7pJOBZyJiZWdvHhHzI6IhIhpqamo6H62VlDOPGcvnZxziUDAbAAoJhBZgTNbrOmBTe30kVQPDgS0dbHsMMEPS42ROQb1N0o+6MX4rA2e8ud6hYDYAFBIIK4AJksZKGkzmIvGSnD5LgDPS8kzg9oiI1D4rVSGNBSYA90TERRFRFxH1aX+3R8R7emE+VqKyQ+Fch4JZUXQaCOmawBxgGZmKoEURsUrSPEkzUrdrgVGSmoHzgblp21XAImA1cAtwbkS80vvTsHLQFgq/ciiYFYUyH+RLQ0NDQzQ1NRV7GNbHfnD341x84yrePmkfrv7nIxlc7b+fNOsJSSsjoqGzfv5JswHnfUfXM6/RRwpm/c2BYAOSQ8Gs/zkQbMByKJj1LweCDWjZoXDO9Q4Fs77kQLABry0UblvjUDDrSw4EKwkOBbO+50CwkuFQMOtbDgQrKQ4Fs77jQLCS876j6/nCrlBY6VAw6yUOBCtJ790VCs84FMx6iQPBSpZDwax3VRd7AGY98d6j6wH4zI2rOOf6lZx93LjiDqifDa4exOS6EVQNyvfoEbOucSBYycsOhdvWPFPcwRTBuJo9+NjUCZx82GgHg/WI73ZqZePRp1+g9cVtxR5Gv3p668t85871rH3qBQeDtavQu506EMxK3M6dwbJVT3Hlrx91MFhevXr7a0nTJK2T1Cxpbp71QyQtTOuXS6rPWndRal8n6cTUNkbSHZLWSFol6WOFT83Msg0aJKa/YT+Wnncc337PkQyuGsTHFtzP2y+/kxvvf5JXdpbOhz4rrk6PECRVAY8AbyfzjOQVwOkRsTqrzznAYRHxYUmzgNMi4t2SJgE3AFOA0cBtwEHA64D9IuJeSa8FVgKnZu8zHx8hmHVu587g1tVPccVtfztiOO9tEzhlso8YKlVvHiFMAZojYn1EbAcWAI05fRqB69LyYmCqJKX2BRGxLSI2AM3AlIjYHBH3AkTEC2QezVlbyMTMrGODBolph776iOHjCzNHDD+/z0cM1r5CAqEW2Jj1uoW//+W9q096BvPzwKhCtk2nl44Alhc+bDPrjIPBuqqQQMh3jJn7ndRenw63lfQa4CfAxyNia943l2ZLapLU1NraWsBwzSxbu8FwmYPBXq2QQGgBxmS9rgM2tddHUjUwHNjS0baSdiMTBtdHxE/be/OImB8RDRHRUFNTU8BwzSyfVwfDUQyudjDYqxUSCCuACZLGShoMzAKW5PRZApyRlmcCt0fmavUSYFaqQhoLTADuSdcXrgXWRMRlvTERMytMJhj2zRsMP7uvhR2v+BYglarTQEjXBOYAy8hc/F0UEaskzZM0I3W7FhglqRk4H5ibtl0FLAJWA7cA50bEK8AxwHuBt0m6P32d1MtzM7MO5AbDkN2q+NeFD/CPl//WwVCh/IdpZga0las+zZW/fpQ1m7cybu89+OjU8Zxy2Giqq3wfzFLmv1Q2s27ZuTP41ZqnueI2B0O5cCCYWY/kBsPYvffgo28bz4zJDoZS40Aws17hYCh9DgQz61VtwXDlbY+y2sFQUhwIZtYnItLFZwdDyXAgmFmfigh+tTpzKsnBMLA5EMysX+QGQ/2o3fno2ybQeLiDYaBwIJhZv3IwDFwOBDMrCgfDwONAMLOicjAMHA4EMxsQIoLb1jzDFbc9wqpNDoZicCCY2YCSLxjmvG0CpzoY+pwDwcwGpNxgOCAdMTgY+o4DwcwGNAdD/3EgmFlJiAh+veYZrvj1Izz8pIOhLzgQzKyk5AuGOW8dz2lH1DoYeqjQQCjoX1nSNEnrJDVLmptn/RBJC9P65ZLqs9ZdlNrXSTqx0H2aWWWRxD9M2odfzDmWa97XwGuHVvPJxQ8y9bI7+XHTRj/BrR90eoQgqQp4BHg70ELmGcunR8TqrD7nAIdFxIclzQJOi4h3S5oE3ABMAUYDtwEHpc063Gc+PkIwqxw+Yug9hR4hVBewrylAc0SsTzteADSSeU5ym0bgc2l5MXCVJKX2BRGxDdiQnrk8JfXrbJ9mVsHajhimHvy6XcHwycUP8pVl6xgxbLdiD6/f/fK8YxlSXdWn71FIINQCG7NetwD/r70+EbFD0vPAqNT+h5xta9NyZ/sEQNJsYDbA/vvvX8BwzaycZAfD7Wuf4cb7N7FjZ+WdPhLq8/coJBDyjSL3PFN7fdprz3e8l/fcVUTMB+ZD5pRR+8M0s3ImiakH78PUg/cp9lDKViEn4lqAMVmv64BN7fWRVA0MB7Z0sG0h+zQzs35USCCsACZIGitpMDALWJLTZwlwRlqeCdwemavVS4BZqQppLDABuKfAfZqZWT/q9JRRuiYwB1gGVAHfjYhVkuYBTRGxBLgW+GG6aLyFzC94Ur9FZC4W7wDOjYhXAPLts/enZ2ZmhfIfppmZlble/cM0MzMrfw4EMzMDHAhmZpY4EMzMDCixi8qSWoEnurn53sCzvTicUuA5V4ZKm3OlzRd6PucDIqKms04lFQg9IampkKvs5cRzrgyVNudKmy/035x9ysjMzAAHgpmZJZUUCPOLPYAi8JwrQ6XNudLmC/0054q5hmBmZh2rpCMEMzPrgAP3LzKQAAAD8klEQVTBzMyAMg0ESd+XtEHS/enr8NQuSd+Q1CzpQUlHFnusPSFpTppLSNo7q73deUo6Q9Kj6euM/HseuCRdL2mdpIclfVfSbqm9bOfcRtJ/Snox6/UQSQvTnJdLqs9ad1FqXyfpxGKMtyckTZV0b/r5vUvS+NReNnNu7+c3a/0bJb0iaWZWW97vZUlHSXoo7e8b6RHGXRcRZfcFfB+Ymaf9JOBmMk9yexOwvNhj7eE8jwDqgceBvTubJ7AXsD79d2RaHlnseXRxzieleQm4AfhIuc85zaMB+CHwYlbbOcC30/IsYGFangQ8AAwBxgKPAVXFnkMX5/sIcHDWPL9fbnNu7+c3rasCbgeWtv0u6+h7mcxzZo5O3/83A9O7M6ayPELoQCPwg8j4AzBC0n7FHlR3RcR9EfF4nlXtzfNE4FcRsSUingN+BUzrvxH3XEQsTfMKMj8EdWlV2c5ZUhXwVeCCnFWNwHVpeTEwNX0ybAQWRMS2iNgANANT+mu8vSSAPdPycP72RMWymXMHP78AHwV+AjyT1Zb3ezl9n+8ZEXenn4sfAKd2Z0zlHAj/kU4dXC5pSGqrBTZm9WlJbeWmvXmWzfzTqaL3ArekpnKe8xxgSURszmnfNbeI2AE8D4yiPOZ8NrBUUguZ/89fTu3lPGcAJNUCpwHfzlnV0fd4S572LivXQLgImAi8kczh1YWpPd95tXKsu21vnuU0/28Cv42I36XXZTlnSaOBfwL+M9/qPG0lP+fkX4GTIqIO+B5wWWov5zm3uQK4MNLTJbP0+dzLMhAiYnM6dbCNzDdT26FjCzAmq2sdfzsULSftzbMs5i/ps0ANcH5Wc7nO+QhgPNAs6XFg9/SoWsiam6RqMqdWtlDic5ZUA0yOiOWpaSHw5rRclnPO0QAsSP+/ZwLflHQqHX+P1+Vp77piX1jpo4s1+6X/ikzafjm9fgevvvB4T7HH2kvzfZxXX1TOO08yR0sbyFyQGpmW9yr2+Ls417OB/wGG5bSX7Zxz5pl9UflcXn2BdVFaPoRXX2BdzwC/wJozx2oyd/Y8KL0+C/hJuc459+c3Z933efVF5bzfy8CK9H3fdlH5pG6Npdj/GL38D7sUGE3m6vxDwMPAj4DXpPUCriZTgfAQ0FDsMfdwnueR+XSwg8wngms6myfwATIX3JqB9xd7Lt2Y8440r/vT18XlPuectuxAGAr8OM3rHmBc1rpPp3+LdXSz4qTI/59PS/8fHwB+0za3cphzZz+/OX13BUJ6nfd7mcxRxcNp/leR7kLR1S/fusLMzIAyvYZgZmZd50AwMzPAgWBmZokDwczMAAeCmZklDgQzMwMcCGZmlvwfiQ6ZK0ziTp0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "tau = 8.\n",
    "\n",
    "def readMatrix(file):\n",
    "    fd = open(file, 'r')\n",
    "    hdr = fd.readline()\n",
    "    rows, cols = [int(s) for s in fd.readline().strip().split()]\n",
    "    tokens = fd.readline().strip().split()\n",
    "    matrix = np.zeros((rows, cols))\n",
    "    Y = []\n",
    "    for i, line in enumerate(fd):\n",
    "        nums = [int(x) for x in line.strip().split()]\n",
    "        Y.append(nums[0])\n",
    "        kv = np.array(nums[1:])\n",
    "        k = np.cumsum(kv[:-1:2])\n",
    "        v = kv[1::2]\n",
    "        matrix[i, k] = v\n",
    "    #化为正负1\n",
    "    category = (np.array(Y) * 2) - 1\n",
    "    return matrix, tokens, category\n",
    "\n",
    "def svm_train(matrix, category):\n",
    "    state = {}\n",
    "    M, N = matrix.shape\n",
    "    #####################\n",
    "    #大于0的化为1\n",
    "    matrix = 1.0 * (matrix > 0)\n",
    "    '''\n",
    "    #构造kernel矩阵\n",
    "    d1 = np.sum(matrix ** 2, axis=1).reshape(-1, 1)\n",
    "    d2 = d1.T\n",
    "    squared = matrix.dot(matrix.T)\n",
    "    dist = d1 + d2 - 2 * squared\n",
    "    k = np.exp(- dist / (2 * (tau ** 2)))\n",
    "    '''\n",
    "    gram = matrix.dot(matrix.T)\n",
    "    squared = np.sum(matrix*matrix, axis=1)\n",
    "    k = np.exp(-(squared.reshape((-1,1)) + squared.reshape((1,-1)) - 2 * gram) / (2 * (tau ** 2)))\n",
    "    \n",
    "    #初始化向量\n",
    "    alpha = np.zeros(M)\n",
    "    #循环次数\n",
    "    n = 40\n",
    "    #系数\n",
    "    L = 1. / (64 * M)\n",
    "    #平均值\n",
    "    alpha_avg = np.zeros(M)\n",
    "    \n",
    "    for j in range(n * M):\n",
    "        #随机取一个样本\n",
    "        i = int(np.random.rand() * M)\n",
    "        #计算函数间隔\n",
    "        margin = category[i] * (k[i, :].dot(alpha))\n",
    "        #grad = M * L * k[:, i] * alpha[i]\n",
    "        grad = L * k.dot(alpha)\n",
    "        if(margin < 1):\n",
    "            grad -= category[i] * k[:, i]\n",
    "        alpha -= grad / ((np.sqrt(j+1)))\n",
    "        alpha_avg += alpha\n",
    "    \n",
    "    alpha_avg /= (n * M)\n",
    "    \n",
    "    state['alpha'] = alpha\n",
    "    state['alpha_avg'] = alpha_avg\n",
    "    state['Xtrain'] = matrix\n",
    "    state['Sqtrain'] = squared\n",
    "    \n",
    "    ####################\n",
    "    return state\n",
    "\n",
    "\n",
    "def svm_test(matrix, state):\n",
    "    M, N = matrix.shape\n",
    "    output = np.zeros(M)\n",
    "    ###################\n",
    "    Xtrain = state['Xtrain']\n",
    "    Sqtrain = state['Sqtrain']\n",
    "    #大于0的化为1\n",
    "    matrix = 1.0 * (matrix > 0)\n",
    "    #做测试集的kernel\n",
    "    gram = matrix.dot(Xtrain.T)\n",
    "    squared = np.sum(matrix * matrix, axis=1)\n",
    "    k = np.exp(-(squared.reshape((-1, 1)) + Sqtrain.reshape((1, -1)) - 2 * gram) / (2 * (tau**2)))\n",
    "    #读取alpha\n",
    "    alpha_avg = state['alpha_avg']\n",
    "    #预测\n",
    "    pred = k.dot(alpha_avg)\n",
    "    output = np.sign(pred)\n",
    "    \n",
    "    ###################\n",
    "    return output\n",
    "\n",
    "def evaluate(output, label):\n",
    "    error = (output != label).sum() * 1. / len(output)\n",
    "    print('Error: %1.4f' % error)\n",
    "    return error\n",
    "\n",
    "def svm(file):\n",
    "    trainMatrix, tokenlist, trainCategory = readMatrix(file)\n",
    "    testMatrix, tokenlist, testCategory = readMatrix('MATRIX.TEST')\n",
    "    \n",
    "    state = svm_train(trainMatrix, trainCategory)\n",
    "    output = svm_test(testMatrix, state)\n",
    "    \n",
    "    return evaluate(output, testCategory)\n",
    "\n",
    "\n",
    "trainMatrix, tokenlist, trainCategory = readMatrix('MATRIX.TRAIN.400')\n",
    "testMatrix, tokenlist, testCategory = readMatrix('MATRIX.TEST')\n",
    "\n",
    "state = svm_train(trainMatrix, trainCategory)\n",
    "output = svm_test(testMatrix, state)\n",
    "\n",
    "evaluate(output, testCategory)\n",
    "\n",
    "size = ['.50','.100','.200','.400','.800','.1400']\n",
    "size1 = [50, 100, 200, 400, 800, 1400]\n",
    "train = \"MATRIX.TRAIN\"\n",
    "error = []\n",
    "for i in size:\n",
    "    file = train+i\n",
    "    error.append(svm(file))\n",
    "    \n",
    "plt.plot(size, error)\n",
    "plt.title(\"error VS szie\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
